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An accurate understanding of the phase diagram of dense hydrogen and helium mixtures is a 
crucial component in the construction of accurate models of Jupiter, Saturn, and Jovian extrasolar 
planets. Though DFT based first principles methods have the potential to provide the accuracy 
and computational efficiency required for this task, recent benchmarking in hydrogen has shown 
that achieving this accuracy requires a judicious choice of functional, and a quantification of the 
errors introduced. In this work, we present a quantum Monte Carlo based benchmarking study of a 
wide range of density functionals for use in hydrogen-helium mixtures at thermodynamic conditions 
relevant for Jovian planets. Not only do we continue our program of benchmarking energetics and 
pressures, but we deploy QMC based force estimators and use them to gain insights into how well 
the local liquid structure is captured by different density functionals. We find that TPSS, BLYP 
and vdW-DF are the most accurate functionals by most metrics, and that the enthalpy, energy, 
and pressure errors are very well behaved as a function of helium concentration. Beyond this, we 
highlight and analyze the major error trends and relative differences exhibited by the major classes 
of functionals, and estimate the magnitudes of these effects when possible. 
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I. INTRODUCTION 

Since the first exoplanet was discovered in 1988, the ad¬ 
vancements in observational techniques and the launch of 
the Kepler telescope have revolutionized our understand¬ 
ing of other solar systems. To date, there are almost 2000 
catalogued extrasolar planetJi^, a large fraction of which 
are Jovian type gas giants and brown dwarves whose 
compositions are more than 90% hydrogen and helium- 
similar to Jupiter and Saturn. Understanding the evo¬ 
lution of these extrasolar planets and their parent star 
systems would thus be greatly aided by accurate mod¬ 
els. Current planetary models for Jovian planets have 
the advantage that by knowing a scant few observables, 
like the mass, radius, luminosity, and element composi¬ 
tions, one can determine the interior structure and time 
evolution of the planet. However, in order to construct 
accurate models, one needs an accurate equation of state 
for hydrogen-helium mixtures at all temperatures, pres¬ 
sures, and species concentrations relevant in planetary 
interiors. 

In the absence of experimental data at these extreme 
conditions, perturbative methodPI, chemical model^and 
ab initio method^^^ have done an excellent job in iden¬ 
tifying phases that could be relevant in planetary mod¬ 
els. For instance, the metallization of dense hydrogen 
liquids is universally believed to occur at temperatures 
and pressures appropriate to the core regions of Jovian 
planets, and should be the source of the planetary dy¬ 
namo responsible for their large magnetic fields. More 
subtly is the possibility of immiscibility of helium in 
dense hydrogen liquids, whereby homogeneous H-|-He 
mixtures driven into the upper atmosphere would con¬ 


dense into helium rich dr oplets , which would rain down 
to the deeper atmospher j^°bi | This process could pro¬ 
vide an additional energy dissipation mechanism as well 
as fundamentally alter the mass distribution of helium 
and heavier elements in the upper atmosphere, and so 
should be treated accurately in planetary models. 

Unfortunately, there is enough quantitative and quali¬ 
tative disagreement among chemical model and ah initio 
based calculations regarding the locations and nature of 
the metallization and immiscibility transitions that plan¬ 
etary scientists routinely treat the equation of state as 
a free parameter in t heir m odels-to be varied to match 
the experimental datcP^^Sl Recent models constructed 
in this manner show an excellent ability to reproduce the 
observed atmospheric depletion of He and excess lumi¬ 
nosity in Saturn, the gravitational m oment s, and esti¬ 
mated ages of both Jupiter and Saturik^li^l, The unfor¬ 
tunate cost of this approach is that it leaves large un¬ 
certainties in other areas of the model-for example the 
core mass and composition, the distribution of heavier 
elements throughout the planet, etc. In principle, these 
uncertainties could be greatly reduced by the additional 
constraints imposed by an accurate equation of state. 

To move forward with the construction of an accurate 
ab initio equation of state with well established error 
estimates, one needs to understand and accommodate 
for two frequently made approximations. The first is 
ideal mixing. The validity of the ideal mixing approx¬ 
imation for the entropy, which has been used in chem¬ 
ical models and ab initio calculations since the 1970’s, 
has only recently been investigated in the context of ab 
initio sim ulat ions using thermodynamic integration (TI) 
techniqued^. There is a current quantitative discrep- 
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ancy of approximately lOOOK in the demixing transition 
at high pressures, and significant qualitative disagree¬ 
ment at lower pressures between the works of Morales 
et al. and Lorenzen et Fortunately, with current 

computational resources, this source of error can be ef¬ 
fectively eliminated through the use of XL 

The second, and least understood source of error is 
in the treatment of electronic correlation effects, typi¬ 
cally through the use of an approximate exchange corre¬ 
lation functional within density functional theory (DFT). 
In pure hydrogen, one can at least quantify the impact 
of density functional errors on the phase diagram, both 
because of explicit benchmarking studies, and because 
the phase diagram has been computed by many differ¬ 
ent groups and methods, ranging from slews of differ¬ 
ent functionals to highly accurate QMC methods. Var¬ 
ious studies have found that the pure hydrogen phase 
diagram is extremely sensitive to the choice of exchange 
correlation functional, primarily because of the presence 
of multiple molecular disassociation and metallization 
phase transitions and crossover JSHIIl contrast, lit¬ 

tle is known of the impact that density functionals er¬ 
rors have on the demixing temperature in dense H-|-He 
mixtures. Since the mid 90’s, PBE has become more 
favored than LDA, but the methodological and quantita¬ 
tive differences between various studies precludes a direct 
comparison. That dense H-l-He mixtures can demon¬ 
strate similar metallization and quasi-molecular transi¬ 
tions should inspire caution. 

Alternative methods for solving the electronic struc¬ 
ture exist. In particular, projector Quantum Monte 
Carlo (p-QMC) is a highly accurate, variational, many- 
body method, well suited to treating electronic corre¬ 
lation in low Z materials. Hydrogen bonding, van der 
Waals forces, and other types of difficult electronic ef¬ 
fects are automatically taken into account. The draw¬ 
back of using QMC in computing entire phase diagrams 
is that it is about two orders of magnitude more expen¬ 
sive than DFT calculations, which makes its widespread 
use in molecular dynamics and ionic Monte Carlo ap¬ 
plications comparatively unattractive with today’s com¬ 
putational resources. Moreover, it is significantly more 
difficult to compute more complex properties relevant to 
planetary physics applications within the QMC frame¬ 
work, like electrical conductivity and viscosity. 

In this work, we instead propose to use projector Quan¬ 
tum Monte Carlo (QMC) methods to benchmark a range 
of density functionals in thermodynamic regimes rele¬ 
vant for helium sedimentation in Jovian planets. Our 
main objective is to identify and understand qualita¬ 
tive error trends and relative differences between vari¬ 
ous classes of density functionals exhibit when used to 
estimate thermodynamic quantities in hydrogen-helium 
mixtures. Though we in many cases strive for quantita¬ 
tive accuracy in addition, this should not distract from 
our main objective. To achieve this, we benchmark the 
errors occurring in the energetics, pressures, and forces 
for each functional, and note how they change as a func¬ 


tion of both density and helium concentration. 

This article is organized as follows. In section |llj we 
first discuss the computational details specific to our cur¬ 
rent work. Then in section [ITT] we present the benchmark¬ 
ing results for global and local energetics, pressures, en¬ 
thalpies, and forces. In section |IV) we explain the error 
trends we observe in terms of the underlying exchange 
functional, after which we conclude. 


II. METHOD 

In this study, we employ the same general methodol¬ 
ogy we used previously in our work on pure hydroger^. 
Thus, we will in this section focus strictly on simula¬ 
tion details and extensions to the method, leaving the 
high level justification and detailed explanation of the 
approach to our previous paper. 


A. Test Sets 


The relevant thermodynamic variables for describing 
the H-l-He phase diagram are the density p, the temper¬ 
ature T, and the helium species fraction xhe-: which is 
defined as: 


Nhc 

- Nh + Nne 


( 1 ) 


In the study of dense hydrogen and helium, it is of¬ 
ten customary to list densities in terms of rg instead of 
p, which are related by H/Ae = |7rrg. is the total 
number of electrons in our system and H is the volume. 

As our goal is to establish the accuracy of density 
functionals around the demixing transition in Saturn 
and Jupiter, all test sets are chosen along the T=7000K 
isotherm. We considered three different densities: rg = 
1.10,1.25,1.34. At each density, in addition to the pure 
hydrogen and pure helium cases, we considered helium 
concentrations between 0 — 20.7%, since these are the 
most relevant compositions for planetary interiors. 

Samples at these thermodynamic conditions consisted 
of charge-neutral cubic cells of 64 electrons with differing 
numbers of H and He ions to ensure charge neutrality. At 
each density and helium concentration, twenty statisti¬ 
cally independent samples were generated from ab-initio 
quantum molecular dynamics simulations in the NVT en¬ 
semble using classical nuclei. The MD simulations were 
performed with the simulation package using 

the PBlP^ exchange correlation functional. 


B. Density Functional Comparison 


For all configurations, we calculated the total energy, 
stresses, and forces using the following functionals in 
Quantum EspresscP^ LDAI^ ( GGA) PBE, revPBE^, 
PBEsoP, BLYPf^, Wu-Coheil^, (metaGGA) M06lPl 




3 


TPS^^ (non-local dispersion corrected) vdW-DFl^, 
vdW-DF^, vdW-DF-C09, vdW-DF2-C0d^, vdW-DF- 
CXP, vdW-optB86]^, and vdW-optB8^. 

For all above functionals, we used a plane wave cut¬ 
off of 800 Ry and a 7x7x7 Monkhorst-Pack grid with an 
offset. We used hard Troullier-Martin pseudopotentialJ^ 
with no core electrons for both H and He that were gener¬ 
ated with OpiunP^ using the PBE functional. To ensure 
no pseudopotential overlap in our test-set, we chose real 
space cutoffs of Xc = 0.37ao and Xc = O.Soq for the H and 
He pseudopotentials respectively. 

We also tested the exact-exchange HSE functional, 
however due to computational and memory limitations, 
we took the following cost saving measures. First, at ev¬ 
ery density and helium concentration we considered, we 
performed only 10 HSE calculations. Secondly, we re¬ 
duced the Monkhorst-Pack grid to 5x5x5 (the same grid 
used to evaluate the Fock operator), and ran the calcu¬ 
lations with VASP because of its compact PAW formal¬ 
ism. A planewave cutoff of 1500eV and 96 bands were 
used for all calculations. We found that this gives the de¬ 
sired accuracy for energy and pressure differences within 
configurations at the same density and helium concen¬ 
tration, however comparisons between different densities 
and helium concentrations might be slightly undercon¬ 
verged. Since we had a limited choice of pseudopoten¬ 
tials, we were unable to perform calculations at = 1.10 
and guarantee the desired level of accuracy for VASP cal¬ 
culations. 


C. Quantum Monte Carlo calculations 

All quantum Mo nte C arlo calculations were done us¬ 
ing the QMCPACKpSlini simulation package. The trial 
wavefunction is taken to be of the single Slater-Jastrow 
form. Singe particle orbitals were obtained from Quan¬ 
tum EspresscP^ using the PBE functional, the same 
Troullier-Martin pseudopotentials described previously, 
and a planewave cutoff of 200Ry. We used pseudopoten¬ 
tials only in the orbital generation step to eliminate the 
electron-ion cusp, which we preferred to handle within 
the electron-ion Jastrow terms. All QMC calculations 
were “all-electron”: we used the bare coulomb interaction 
between electrons and electrons, and between electrons 
and nuclei. 

Eor the Jastrow factor, we used short-ranged one-body 
and two-body functions of b-spline form. Eor H and 
He, the one-body terms were spin-independent. The one 
body term for each species was a sum of two functions. 
The first was a “core” jastrow, which had a real space cut¬ 
off of Xc = l.Ooo, 8 knots, and had the suitable electron- 
ion cusp condition imposed. The second had a cutoff of 
L/2 with 8 knots and no cusp-condition imposed. For 
the two-body functions, we separately included same- 
spin and opposite-spin e-e terms, each with a cutoff of 
Tc = T/2 and correct cusp conditions imposed. 

Our wavefunctions were optimized with the linear 


method^. After obtaining a good initial guess for the 
jastrow parameters from a single r^ = 1.10 configuration 
with 4 He atoms, all variational parameters were simulta¬ 
neously optimized using an initial variance minimization 
step, followed by 10 energy minimization steps. Conver¬ 
gence of the minimization procedure was checked. 

Energies, pressures, and the structure fact or wer e cal¬ 
culated using Reptation Monte Carlo (RMC)PIIS1_ Qur 
target statistical error bars for the energies and pres¬ 
sures were 0.008 mHa/electron and 0.3GPa respectively. 
For all but the pure helium configurations, we used a 
time-step of r = 0.0075iJa“^ and projection time of 
(3 = 4.5i7a“^. These choices were found to yield time- 
step and mixed-estimator errors for the potential en¬ 
ergy that were comparable to the desired error bars. 
For the pure helium configurations, we fixed the projec¬ 
tion time at /3 = 4.5i7a“^ and ran with time steps of 
T = 0.0075i7a“^ and r = 0.0037577a“^. All accumu¬ 
lated quantities were then linearly extrapolated to zero 
time step. 

Forces were computed using the Chiesa, Ceperley, 
Zhang estimator^ adapted to periodic boundary condi¬ 
tions, which we detail in the supplemental information. 
We used a real-space cutoff of 7^ = l.Ooo and a smooth¬ 
ing polynomial of degree M = 3. Based on several statis¬ 
tical tests detailed in the supplemental information, we 
found that this choice of parameters yielded systematic 
errors that were less than the error bar on the hydrogen 
force components, which was approximately 2mHa/bohr. 
This resolution was sufficient to clearly distinguish differ¬ 
ent functionals. We used diffusion Monte Carlo with a 
time-step of r = 0.01i7a“l and a population size of 512 
walkers, which we found converged the local energy (but 
not the potential) to within error bars. To correct for 
the mixed-estimator problem we used extrapolated force 
estimates. All systematic errors are expected to be less 
than the statistical error bar. 

We applied the following finite size corrections. To 
reduce shell effects, we used canonical twist-averaged 
boundary conditions (TABC) on a 4x4x4 Monkhorst- 
Pack gricP^. For the potential energy correction, we used 
the leading order Chiesa correction based on pure esti¬ 
mates of See(k^^- For the kinetic energy correction, we 
detail this in the supplemental information. Since most 
configurations are in the metallic state, there is also a 
kinetic energy error arising from the fact that we are 
attempting to reproduce a fermi-surface with only 64 
electrons. To correct for this, we used the PBE func¬ 
tional and estimated the energy error between a twist- 
averaged unit cell and a 7x7x7 MP grid. This correction 
scheme was tested against several supercell calculations 
at Xs = 1.10 and Xg = 1.34 across all helium concen¬ 
trations. We expect the absolute energy errors (across 
all densities and helium concentrations) to be less than 
0.5mHa/electron, and the pressure errors to be approxi¬ 
mately IGPa. Details are in the supplemental informa¬ 
tion. 

Regarding QMC force errors, a few caveats are worth 
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mentioning. No finite-size correction scheme beyond 
twist-averaging was applied to the forces. Though this 
should eliminate the bulk of the finite-size error, there is 
still a residual finite size error coming from the fact that 
p(r) is not quite converged to the thermodynamic limit. 
By comparing the forces from a twist-averaged KZK cal¬ 
culation in the unit cell with those of a close to converged 
7x7x7 MP grid, we estimate that the standard deviation 
of the residual finite-size error in the force estimates is 
approximately 2mHa/bohr. 

D. Error Analysis 

1. Scalar Quantities 

As in our previous paper, we consider a test set S 
with M configurations, {Rq ... Rm}- For each configu¬ 
ration Ri, we define a density functional error 5vA(Ri) = 
where A is some observable (e.g. 
total energy, pressure), “DF” is the density functional, 
“QMC” is the QMC reference value. 

In addition to defining average errors over a test set S, 
which we denote {6A^^)s, we define a general class of 
shifted mean absolute errors as: 

= E (2) 

^ RiGS 

Here, is an density functional dependent offset. The 
standard “mean absolute error” corresponds to the choice 
of = 0, which we will label as (|(5A|)s. However, 
we will in this paper find cause to use other choices of 
for different observables, specifically for measures of 
“global” and “local” energetics, which we will take care 
to explain as they arise. 

2. Forces 

Let ii denote the force on ion i, and F = {fi ,f2, . . . , fAf} 
is the 3N dimensional vector of all ionic force compo¬ 
nents. Because of the large number of force components, 
we spend this section describing how we reduce the di¬ 
mensionality of this data set to construct a handful of 
force-error measures. 

One of the simplest measures of force errors we can 
devise is the mean absolute force error (|5f^^|)s. In¬ 
tuitively, this is the ensemble average magnitude of the 
force error vector Sfi = 

We can do significantly better than establishing which 
functional has better forces on average. Indeed, know¬ 
ing the force on each atom gives significant insight into 
the system’s local structure. Rigorously, one can do this 
through calculating the “potential of mean force” w{r), 
which is directly related to the pair correlation function 
by g{f) = e'xjp{—j3w{r)). Where the proper distribution 
is not known or hasn’t been sampled adequately, one can 


often approximate this quantity fairly well and repro¬ 
duce the major features of the pair correlation function, 
as is routinely done in force-matching with classical po¬ 
tentials. 

For each density functional, we will define the following 
error measure relative to QMC. Let a and ji, denote 
two particles of species pt and v respectively. Denoting 
“ ’’v ~ define , 

( 3 ) 

Based on this definition, if is positive (nega¬ 

tive), it overbinds (underbinds) species of type /i and v 
at a distance r. 

Note that we use in this definition, since our 

configurations are sampled from QMD using the PBE 
functional. However, if we could replace with 

EQMC^ then {Sff^J^^ir))QMC would be related to the 
density functional error in the potential of mean force 
(r) (relative to the QMC distribution) by 

{Sfl^lAr))QMC ={r) (4) 

In any case, given that E^^^ and E^^^ produce 
qualitatively similar distributions of ionic configurations, 
'^bl enable us to see over which regimes a 
density functional overbinds or underbinds, which would 
give strong indications as to how the g{r) would change 
based on density functional. 

III. RESULTS 
A. Global Energetics 

Suppose we are interested in assessing to what extent 
the average error in the total energy changes as a function 
of helium concentration. Depending on the error scaling, 
this can affect the both the Helmholtz and Gibbs free- 
energies of mixing, and thus the location of the H-|-He 
immiscibility transition. 

To measure this quantity, we define a measure of 
“global energetics” as follows: for a given p, we build 
an aggregated test set S'{p) which consists of all test 
sets at all helium concentrations with a given electronic 
density p. We then choose c^^{p) to be the median of 
{6E^^}s'{p)- With this definition for we define 

the “global energetic error” of the test set S{p,XHe) by 
using Eq. which for shorthand we will refer to as 

{\s^n)9,s- 

In Fig. we show the global energetic error for all 
functionals at three densities, averaged over all helium 
concentrations. We see that the best performing func¬ 
tionals over all densities are the meta-GGA functionals 
TPSS and M06-L respectively, with global energetic er¬ 
rors that are under half that of PBE. After these, the 
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FIG. 1: {\5E^^\)global averaged over all helium concentra¬ 

tions for all considered functionals. The different bar col¬ 
ors/patterns denote the different densities. 


best performing functionals are the semi-local GGA’s 
BLYP and revPBE, followed by optB86b-vdW and vdW- 
DF. The worst performing functionals are LDA, HSE, 
PBEsol, and WC, with global energetic errors approxi¬ 
mately twice that of PBE. Though PBE has better than 
average performance in this regime, one can gain a de¬ 
cent fraction of a milihartree in accuracy by switching to 
a metaGGA or a properly tuned GGA. 



Xjjo (%) 


FIG. 2: {5[E{xHe) — E{xHe = 1)]) vs. XHe for all considered 
functionals at (left) Vs = 1.10, (middle) Vs = 1.25, and (right) 
Vs = 1.34. All energies are measured relative to the aver¬ 
age energy of all pure helium conhgurations at the specihed 
density. 

It is also useful to consider how accurately the energy 
difference between systems with different helium con¬ 


centrations are captured with a density functional. For 
specificity, we look at {S[E(xHe, p) — E[xHe = IjP)]) in 
Fig.i for all considered functionals and helium concen¬ 
trations. What we see is that relative to the pure helium 
configurations, BLYP, vdW-DF, and vdW-DF2 overes¬ 
timate the energy difference between the mixed hydro¬ 
gen/helium configurations, whereas all other function¬ 
als underestimate this difference. Additionally, though 
all curves exhibit noticeable nonlinearity in the xue = 
0 — 20% range, almost all curves have magnitudes that 
monotonically decrease to zero. The exception is M06L, 
which appears as though it reaches a maximum energy 
error somewhere between xrb = 20 — 100%. 


B. Local Energetics 

Even if the total energy error for an DFT-MD simula¬ 
tion at specific density and helium concentration averages 
to a small value based on the previous error measure, it is 
still possible for energy differences between similar con¬ 
figurations to have large errors. To measure the spread of 
the error distribution, we use the following. For each test 
set S{p,XHe) corresponding to a given density p and he¬ 
lium concentration xhb, we set to be the median of 
the set {SE°^}s(p^ooh.)- Using Eq. I with this choice of 
test set and gives our “local energy” measure, which 
we denote as {\6E^^\)i^s- Note that this is similar to the 
“global energetic error” measure used previously, except 
now the test set S and the reference set S' are the same. 



FIG. 3: {\SE^^\)local averaged over all helium concentra¬ 

tions for all considered functionals. The different bar col¬ 
ors/patterns denote the different densities. 

In Fig. we show {\SE^^\)iocai at three different 
densities for all functionals considered. We note that the 
characteristic error scale is now 0.25 mHa/electron, and 
that the differences between density functionals is signif¬ 
icantly less pronounced than in the global energetic case. 
So small are the differences, in fact, that little can be 
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said for or against the vast majority of functionals by 
this metric, but the outliers are worth mentioning. LDA, 
PBEsol, vdW-DF2, and M06L are the worst performers, 
whereas HSE is perhaps the best. These trends are con¬ 
sistent with the global energetic trends, with the notable 
exception of M06L, which accurately captures global en¬ 
ergetics, and HSE, which exhibits poor global energetic 
performance. 


C. Pressures 

For all test sets and functionals, we computed {SP^^)s 
and {SP^^)s- As observed in pure hydrogen, we found 
the mean pressure errors to be systematically off, whereas 
over all helium concentrations and densities is 
virtually indistinguishable from the statistical error bars. 
{SP^^)s on the other hand can be quite sizable. For the 
rest of this section, we will plot relative mean pressure 
errors instead of absolute errors because of the dramatic 
change in pressure as one changes the helium concentra¬ 
tion. For example, the pressure drops at Vg = 1.10 from 
over IT Pa with pure hydrogen to around SOOGPa for 
pure helium. 

In Fig. we plot {5P^^)s/{P^^^) (averaged over 
all helium concentrations) at three different densities for 
all functionals considered. The trend observed is very 
much the same one observed in our previous benchmark¬ 
ing studies of pure hydrogen: accurate energetics are 
compensated by poor pressure estimation. For exam¬ 
ple, PBEsol, Wu-Cohen, LDA, though the worst per¬ 
formers energetically, provide the most accurate pressure 
estimates, missing the correct pressure estimate by less 
than 1%. The worst shown functionals happen to be 
vdW-DF and BLYP, which were known for their accu¬ 
rate energetics. vdW-DF2 is a bit of an exception, in 
that it has poor energy and poor pressure estimation. 
Note that TPSS and M06-L functionals are not plotted. 
This is because relative to all other functionals, the pres¬ 
sure errors are quite significant. Averaged over all three 
densities, M06-L and TPSS have average pressure errors 
of approximately -31% and -17% respectively, both ex¬ 
hibiting increasingly poor performance as the density is 
decreased. 

In Fig. we plot the mean pressure error versus he¬ 
lium concentration for all considered functionals. Gener- 
ically, though some functionals might underestimate or 
overestimate the pressure errors as the helium fraction 
is increased, the dependence of the pressure errors on 
XHe is very well behaved. The magnitudes for almost 
all functionals increases monotonically as Xhc increases, 
reaching its maximum error for pure He. The exceptions 
are OLYP, where the pressure error has a positive slope 
and changes sign at some nonzero xhe, and WC, which 
reaches a maximum at some nonzero XHe and then de¬ 
creases towards a minimum at XHe = 1- Though we did 
not investigate helium concentrations higher than about 
20%, the smoothness of the pressure errors as a func¬ 


tion of XHe is reassuring, as it opens up the possibility 
of fitting these errors and correcting for them in post 
processing. 



FIG. 4: in units of (%) averaged over all 

helium concentrations for all considered functionals. The dif¬ 
ferent bar colors/patterns denote the different densities. Not 
shown: M06-L and TPSS. 



FIG. 5: {6P°^)/{P'^^^') in units of (%) vs. XHe for all 

considered functionals at (left) = 1.10, (middle) Vs = 1.25, 
and (right) Vs = 1.34. The different colors/shapes are given 
in the legend and denote the density functional. Note the 
broken axis, which shows the mean relative pressure error for 
the pure helium configurations. 


D. Enthalpies 

The mean enthalpy error g is given by 

{6H^^)s = {5E^^)s + V{6P^^)s, which we can com- 
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bine with the results from |III B| and |III C to study the 
errors in the predicted DFT enthalpy. 

In Fig. 1^ we show the the average enthalpy er¬ 
ror versus helium concentration for = 

1.10,1.25,1.34. At fixed density p, is measured rel¬ 
ative to the enthalpy of pure helium at density p. What 
we see is for the most part qualitatively similar to what 
we saw for the energy errors in Fig. However, we see 
that the tested functionals possess varying degrees of er¬ 
ror cancellation. Some functionals definitely benefit from 
error cancellation: specifically vdW-DF, LDA, and HSE. 
This can reduce the absolute enthalpy by as much as 3-4 
mHa/electron depending on the functional and density. 
Others suffer from error addition, such as PBE and most 
dramatically OLYP, which has an enthalpy error almost 
lOmHa/electron higher than its corresponding energy er¬ 
ror. Lastly, there are some functionals which exhibit nei¬ 
ther error cancellation nor addition, namely BLYP and 
all the newer vdW functionals. 



FIG. 6: {S[H^^{xHe) - = 1)]) vs xue for all 

considered functionals at (left) = 1.10, (middle) = 1.25, 
and (bottom) Xg = 1.34. All energies are measured relative 
to the average energy of all pure helium configurations at the 
specified density. Note that the reference enthalpy for each 
density is taken to be the mean enthalpy of the pure helium 
configurations at that density. 


E. Forces 


1. Total Force Errors 

The most natural question we can ask ourselves is 
which functional has the most accurate forces on aver¬ 
age. To this end, in Fig. we compute (|5f^^|) over 
all atoms and helium concentrations. Before interpreting 



FIG. 7: {\5t[j^\) aggregated over all helium concentrations 

for all considered functionals. The different colors denote dif¬ 
ferent densities. 


Fig. 0 however, we mention an important caveat. Be¬ 
cause of the large statistical noise relative to the individ¬ 
ual force components, any mean absolute error measure 
is going to be saturated by the statistical noise. Thus, 
the absolute magnitude of the mean absolute force er¬ 
ror should not be interpreted as being indicative of the 
underlying functional performance. In spite of this how¬ 
ever, the fact that the statistical noise introduced in each 
force error measure is identical across all functionals, we 
are able to compare different functionals relatively and 
establish which ones are better and by how much. 

With the above caveats in mind, we note in Fig. 0 
that we observe the trend in mean absolute force errors 
is extremely similar to the trends we saw for local ener¬ 
getic errors. Specifically, HSE is among the best, whereas 
vdW-DF2, LDA, and M06-L are among the worst. There 
are some slight differences however. TPSS and optB86b- 
vdW seem to perform better by force error measures than 
the local energetics would suggest. Additionally, though 
BLYP and vdW-DF still seem to outperform PBE, the 
difference is not nearly as apparent functionals as the 
local energetics would suggest. 


2. Local Force Errors 


In this section, we tabulate the average density func¬ 
tional force errors as a function of distance (^/,f^(r)), as 
described in Section [IlD 2| We consider H-H, H-He, and 
He-He interatomic forces. In principle, these force errors 
will not just depend on density functional, but also on 
density and helium concentration. We will address these 
later two points first, as they will greatly simplify the 
analysis. 

The first question is how does the average force change 
as a function of density. We show in Eig. I 
versus r/r^ calculated using the PBE functional. We 



































































consider H-H (top), H-He (middle), and He-He (bottom) 
forces. The first two were calculated at a helium fraction 
of 20.75%, whereas the last was calculated in pure helium 
for statistical reasons. For each plot, we overlay the plots 
of at the densities = 1.10,1.25,1.34. What 

is obvious is that with the exception of the sometimes 
large statistical fluctuations, no significant quantitative 
or qualitative differences exist in the mean force errors 
versus distance. 


The second question is how does the average force 
change as a function of helium concentration. In Fig. 
^ we show the same general plots of using 

me PBE functional as in Fig. but this time overlaying 
plots of different helium concentrations instead of dif¬ 
ferent densities. All plots were calculated at a density of 
Tg = 1.25. Note that within error bars, {5f^^^{r)) shows 
a remarkable insensitivity to helium concentration. 


Given the insensitivity of the average force errors for 
PBE to both density and helium concentration, we plot 
considered functionals. The helium 
fraction was chosen to be 1.6% for the H-H (top), 20.75% 
H-He (middle) plots, and 100% for the He-He plot. The 
dens ity wa s chosen to be Vg = 1.25. Recalling from Sec¬ 
tion HD2 that > 0 implies overbinding rela¬ 

tive to GMC. 


For the H-H forces at the top of Fig. the BLYP, 
vdW-DF, and vdW-DF2 functionals all exhibit a strong 
propensity to overbind in the 1 < r/vg < 1.5, with vdW- 
DF2 overbinding the most. TPSS overbinds the least in 
the region 1 < r/rg < 1.2 but then under binds slightly 
up to r/vg = 2.2. All other functionals underbind in the 
region 1 < r/vg < 1.5, with HSE underbinding the least 
and LDA the most. Though its hard to tell with the 
noise, HSE has the lowest absolute error in the region 
1 < r/rg < 1.5, followed by optB86b-vdW, vdW-DF 
and BLYP, and then by the combination revPBE, PBE, 
vdW-DF-CX, vdW-DE-C09, and vdW-DF2-C09. 


For the H-H forces at the top of Fig. there seem to 
be three distinct regions in space, whose boundaries are 
roughly where the vast majority of DF errors cross the 
r-axis. I will refer to these as region I (1 < r/rg < 1.5), 
region H (1.5 < r/rg < 2.2), and region HI {r/rg > 2.2). 
These should roughly correspond to the first, second, and 
third coordination shells. Notice that with the excep¬ 
tion of TPSS and M06-L, if a functional overbinds in 
region I, it will almost certainly underbind in region H, 
and overbind again in region HI. This is not entirely un¬ 
expected. The ion-ion force depends only on electron 
densities, and so if two protons overbind because of an 
increased electronic charge between them, this decreases 
the electronic charge elsewhere, leading to underbinding 
in the charge depleted region. 

There are only a few functionals that overbind the H- 
H interaction in region I: BLYP, vdW-DF, and vdW- 
DF2, and TPSS (only ior r/rg « 1). The rest underbind, 
though to varying degrees. If we try to determine which 
functionals have the smallest error magnitudes in region 
I, we find that the trend is very similar to what we saw 


before in the mean absolute force and local energetic sec¬ 
tions. HSE, vdW-DF, BLYP, TPSS, and optB86b-vdW 
have the smallest errors in regions I, though further dis¬ 
crimination is difficult given the error bars. In region 
H on the other hand, HSE and optB86b-vdW seem to 
have measureably smaller error magnitudes than vdW- 
DF, BLYP, and TPSS. 

For the H-He forces in the middle of Fig. the dif¬ 
ferences between different functionals are more striking. 
HSE and TPSS have the best average performance in 
the region 1.5 < r/rg < 2.0. However, BLYP also per¬ 
forms exceptionally well, slightly underbinding hydrogen- 
helium pairs by less than ImHa/bohr. The worst per¬ 
forming functionals are LDA, which overbinds the H-He 
interaction, and vdW-DF2, which underbinds. 

Lastly, we consider the He-He forces at the bottom of 
Eig. |10[ the error bars are somewhat large, but there are 
some obvious trends still visible. All functionals overbind 
the He atoms, although LDA overbinds the most. The 
functionals that underbind the least are either vdW-DF2, 
HSE, or TPSS, followed by vdW-DF and then maybe 
BLYP. 


IV. DISCUSSION 
A. Energeties 

In this section, we will try to reconcile the differ¬ 
ences we observed between the local and global energetic 
trends. We believe that the trends observed in the lo¬ 
cal energetic errors are mostly described by the impact 
that the density functional has on the charge density and 
forces, we will spend more time discussing that point in 
section |IVD[ Por now, we will try to grapple with how 
little bearing the local energetic errors had on the global 
energetic errors. 

The main point to realize is that the global energetic 
errors are going to be dominated by errors in the en¬ 
ergy differences between configurations with different he¬ 
lium concentrations. Calculating these types of energy 
differences accurately means one of two things: having 
small total energy errors, having error cancellation, or 
both. Error cancellation is possible only if the H-H, H- 
He, and He-He interactions are described in roughly the 
same way. But from our force discussion, we saw that 
many functionals will overbind one type of interaction 
while underbinding another. An instructive compari¬ 
son would be between vdW-DF and BLYP. Both have 
very similar performances for forces and local energet¬ 
ics, with vdW-DF having a slight edge on both. How¬ 
ever, BLYP has noticeably smaller global energetic er¬ 
rors than vdW-DP. Consider, for simplicity, the error in 
E{xHe = 0) ~ E{xHe = !)■ Looking at the force er¬ 
rors, BLYP and vdW-DE overbind the H-H interaction 
in almost the exact same way. However, vdW-DF under¬ 
binds the He-He interaction noticeably less than BLYP 
does. The consequence is that the energy of the sampled 
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pure He configurations (from PBE based QMD) are go¬ 
ing to be higher than vdW-DF will predict. However, 
both vdW-DF and BLYP are overestimating the energy 
of a pure sampled H configuration by about the same 
amount, so the magnitude of 5[E{xHe = 0) — E{xHe = 1)] 
will be greater for vdW-DF than it will be for BLYP be¬ 
cause of error cancellation. This argument only goes so 
far however, as there are other contributions to the en¬ 
ergy than just the electrostatic contribution. HSE for 
instance, though probably describing the charge density 
very well and having small local energetic errors, has very 
poor global energetics. 


B. Pressures 

As in our previous hydrogen benchmarking work, we 
observed that favorable energy errors were anti correlated 
with favorable pressure errors-the most dramatic cases 
being LDA and vdW-DF/BLYP. Though this might seem 
paradoxical, the recognition of a tradeoff between accu¬ 
rate bulk moduli and lattice constants versus cohesive 
and atomization energies is well known. In fact, Perdew 
et al. argue that this tradeoff is a necessary consequence 
of the limited form of the GGA functionaP^. Basically, 
to reproduce accurate binding and atomization energies 
necessitates a fj, that is higher than that associated with 
the gradient approximation expansion, which is neces¬ 
sary to recover the slowly-varying electron gas limit. We 
find that almost all the best performing functionals in 
this regard are still the functionals that have been de¬ 
signed with these constraints in mind: specifically, the 
LDA, Wu-Gohen, and PBEsol functionals. 

The hybrid functional HSE is the exception, instead us¬ 
ing the reduced self-interaction error to achieve a better 
estimate of the pressure. Not only did we observe rea¬ 
sonably accurate pressure estimation, HSE additionally 
had some of the smallest errors with local energetics and 
forces. The global energetics errors were among the worst 
tested, but this might be due more to a more inequitable 
treatment of H and He and a lack of error cancellation, 
rather than a consequence of large absolute errors. 


C. Enthalpies 

When constructing the equation of state for H-|-He 
mixtures, the most important thing we have to worry 
about is having accurate enthalpies. Though the entropy 
term is also important, is far less sensitive to the choice of 
density functional than the energies and pressures. That 
being said, one can cut the enthalpy errors by approxi¬ 
mately 50-60% relative to PBE (from llmHa/electron to 
4mHa/electron in pure hydrogen) by using either BLYP 
or vdW-DF. Improving the enthalpy errors beyond this 
without using some sort of post-processing scheme might 
be somewhat difficult. The 4mHa for vdW-DF and BLYP 
is in large part due to signihcant (though noticeably in¬ 


complete) error cancellation. Given the inherent tradeoff 
between energy errors and pressure errors discussed pre¬ 
viously, one should be extremely careful correcting each 
piece individually, especially if one can’t fall back on a 
higher level of theory to verify. 

D. Forces 

From our analysis of local energetic errors and forces, 
we saw that there is a strong though not perfect correla¬ 
tion between small energetic errors and small force errors. 
As mentioned before, having accurate electron-ion forces 
depends only on the ability to accurately reproduce the 
electronic charge density, whereas the local energetic er¬ 
ror measures additionally fold in errors in the treatment 
of electron-electron correlation effects. 

With this in mind, the superior performance of the 
HSE functional in minimizing the local energetic er¬ 
rors and force errors most likely stems from its ability 
to produce a reasonable charge-density. This shouldn’t 
be surprising, as the introduction of exact exchange fa¬ 
vors charge localization through the reduction of self¬ 
interaction errors. After HSE, TPSS seems to produce 
reasonable charge densities. Among the GGA’s and vdW 
corrected GGA’s, the vdW-DF, BLYP, and optB86b- 
vdW functionals seem to produce reasonable charge den¬ 
sities, most likely because the underlying exchange func¬ 
tionals are skewed to energetically favor bonding and 
charge localization. 

E. Role of Exchange 

After benchmarking several GGA, meta-GGA, hybrid, 
and non-local vdW functionals, it is perhaps safe to say 
that most of the differentiation in performance between 
functionals we considered stems from the treatment of 
exchange, and not from the addition of sophisticated non¬ 
local correlation effects. This conclusion follows from two 
pieces of evidence. The first is that as far as global and 
local energetic errors are concerned, the two best per¬ 
forming functionals are vdW-DF, a non-local vdW func¬ 
tional, and BLYP, a GGA. Beyond just having compara¬ 
ble performance, the total magnitude and scaling of local 
and global energy errors with helium concentration are 
very similar. One would expect that if vdW type correla¬ 
tion were necessary for an accurate description of dense 
H-|-He, that there wouldn’t be any GGA’s performing 
nearly as well as the non-local vdW functionals. 

The other indication that the improved energetics is 
driven by the exchange piece comes from comparing the 
performance of the non-local van der Waals function¬ 
als. vdW-DF, vdW-DF-C09, and vdW-DF-GX all use 
the same non local correlation functional, differing in 
their choice of exchange functional only. The same is 
true of vdW-DF2 and vdW-DF2-G09. We found that 
vdW-DF-G09 and vdW-DF2-G09 were virtually indistin- 
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guishable energetically, indicating the small role played 
by the difference in the van der Waals correlation piece. 
However, there is a significant difference between vdW- 
DF-C09 and vdW-DF, or between vdW-DF2-C09 and 
vdW-DF2, each pair demonstrating either a propensity 
underbind or overbind respectively relative to QMC. 

It turns out that the best performing density func¬ 
tionals exhibit some common trends in their underly¬ 
ing exchange functionals. The exchange functionals for 
GGA’s are given by Ex[p\ = f drp(r)e^°”^(r)Fx(s(r)), 
where is the Slater-type exchange for the homoge¬ 
neous electron gas, is the “enhancement factor”, and 
s = |Vp|/[2(37r^)^/^p^/^] is the “reduced density gradi¬ 
ent”. Before getting into similarities in F^ responsible 
for decent or poor energetic or pressure performance, we 
need to know which values of s are relevant in our system. 
After analyzing the PBE and BLYP charge densities for 
a single sample configuration from each density and he¬ 
lium concentration, we conclude that s is bounded by 
0 < s < 1.8 for all configurations of interest (s < 0.8 for 
pure H configurations). Unsurprisingly, the largest gra¬ 
dients occur in pure helium configurations at low density, 
whereas the smallest gradients occur in pure hydrogen at 
high density. 

Within the semilocal GGA functionals, we can ex¬ 
plain better or worse energetic performance relative to 
PBE by looking at F^. We saw that the global ener¬ 
getic, local energetic, and force errors followed the pro¬ 
gression of decreasing accuracy: BLYP, revPBE, PBE, 
and PBEsol. Looking at the underlying enhancement 
factors (BLYP uses B88 exchange), we see the following 
trend: Ff^^ > > F^^^ > Ff for all “s” 

in the relevant range for hydrogen. prevPBE ^ 
from about s = 0.8 onwards (they cross again at much 
larger s), but this doesn’t affect the description of hy¬ 
drogen. This implies that the best performing func¬ 
tionals for energies and forces are working by lowering 
the energy contributions coming from the larger reduced 
gradients, which favors charge localization and bond¬ 
ing. Additionally, noting that the relative difference 
Ef®® - « 0.005 at s = 0.4, and recalling how 

much the H-H forces and local energetics changed with 
respect to functional implies that the electronic struc¬ 
ture around protons is very sensitive to the treatment 
of exchange at these densities. Looking at how similar 
the local energetic errors and He-He forces were for pure 
helium configurations for different functionals would in¬ 
dicate that the helium is not nearly as sensitive to the 
choice of exchange functional. 

One can perform the same type of analysis with the 
vdW-DF type functionals. vdW-DF, optB86b-vdW, and 
vdW-DF-GX use the revPBE, optB86b, GX exchange 
functionals respectively. We previously saw that for ener¬ 
getic and force errors, the progression towards decreasing 
accuracy follows the sequence vdW-DF, optB86b-vdW, 
and vdW-DF-GX. Looking at the underlying enhance¬ 
ment factors, we find that > E°ES866 > pCX ^ 

vdW-DF and optB86b-vdW perform comparably, but 


vdW-DF overbinds relative to QMG whereas optB86- 
vdW underbinds. We forgo a direct exhange functional 
comparison between the vdW functionals and the GGA’s, 
primarily because of the “exchange consistency” compli¬ 
cation stemming from the use of a different “outer” and 
“inner” exchange correlation functional. 

Deeper relationships between the functional form of 
Fx{s) and corresponding errors can be deduced from the 
previous discussion. However, we leave these consider¬ 
ations to future publications, since our current focus is 
on hydrogen-helium thermodynamics and not on density 
functional development. 


V. CONCLUSIONS 

In this paper, we have used projector Quantum Monte 
Carlo to benchmark a collection of the most popular den¬ 
sity functionals, ranging from GGA, to non-local disper¬ 
sion corrected, to meta-GGA. We were able to quantify 
the errors for most quantities that are relevant for con¬ 
structing an equation of state: specifically the pressures, 
local and global energy differences. As a result of our 
analysis, we can conclude that significant reduction of en¬ 
thalpy errors and a much better description of hydrogen 
helium interactions can be attained by using the TPSS 
metaGGA, the BLYP GGA, or the nonlocal vdW-DE 
functionals. 

Beyond just identifying the most accurate density func¬ 
tional and quantifying its errors, we have demonstrated 
the common features of the best performing functionals, 
specifically in the shape and limiting behavior of the en¬ 
hancement factors for the exchange functionals. The un¬ 
derlying exchange pieces for both vdW-DF, BLYP, and 
revPBE tend to emphasize bonding in the energetics, 
which is well known in the DFT literature. The impor¬ 
tance of this work is that it specifies quantitatively just 
how important this bonding character is for an accurate 
description of dense hydrogen helium mixtures. Know¬ 
ing this, and how much the various exchange correlation 
functionals overbind or underbind, should facilitate the 
optimization and deployment of new functionals for map¬ 
ping out the H-|-He phase diagram. 
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FIG. 8: {Sf^-^{r)) vs. r/va as density is changed. The dif¬ 

ferent marker colors/styles represent different densities. (Top) 
{Sfu-Hi''')) calculated at life = 20.7%, (middle) {Sfu-uAr)) 
calcnlated at xh& = 20.7%, (middle) {Sfue-Htti r)) calculated 
at xh£ = 100% 
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FIG. 9: {5f^F^{r)) vs. r/r^ as helium concentration is 

changed. The different marker colors/styles represent dif¬ 
ferent helium concentrations. (Top) (middle) 

('5/if-fe(»')), (middle) {S/nf-Heir))- All configurations are 
at a density of Vs = 1.25. 
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FIG. 10: {5f^f^{r)) vs. r/vs as the functional is changed. 

The different marker colors/styles represent different den¬ 
sity functionals. (Top) {5-%{''')) a-t = 1.6%, (mid¬ 
dle) (SfH^Heir)} at XHe = 20.7%, (middle) {Sfue-Heir)} at 
XHe = 100%. All configurations are at a density of Vs = 1.25. 




























16 



FIG. 11: Plot of enhancement factors Fx{s) for all the best 
performing functionals. 







